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We consider the Bak-Tang-Wiesenfeld sandpile model on square lattices in different dimensions 
(_D < 6). A finite size scaling analysis of the avalanche probability distributions yields the values of 
the distribution exponents, the dynamical exponent, and the dimension of the avalanches. Above 
the upper critical dimension Du ~ 4 the exponents equal the known mean field values. An analysis 
of the area probability distributions indicates that the avalanches are fractal above the critical 
dimension. 

PACS number: 05.40. -fj 



I. INTRODUCTION 

Bak, Tang and Wiesenfeld introduced the concept 
of self-organized criticality (SOC) and reahzed it with 
the so-called 'sandpile model' (BTW model). The steady 
state dynamics of the system is characterized by the prob- 
ability distributions for the occurrence of relaxation clus- 
ters of a certain size, area, duration, etc. In the crit- 
ical steady state these probability distributions exhibit 
power-law behavior. Much work has been done in the 
two dimensional case. Dhar introduced the concept of 
'Abelian sandpile models' which allows to calculate the 
static properties of the model exactly ||^, e.g. the height 
probabilities, height correlations, number of steady state 
configurations, etc (^^. Recently, the exponents of the 
probability distribution which describes the dynamical 
properties of the system were determined numerically 
|3|. On the other hand both mean field solutions (see 
|7| and references therein) and the solution on the Bethe 
lattice 1^ are well established and both yield identical 
values of the exponents. The mean field approaches are 
based on the assumption that above the upper critical 
dimension the avalanches do not form loops and the 
avalanches propagation can be described as a branching 
process § . Despite various theoretical and numerical ef- 
forts the value of D„ is still controversial. In an early 
work, Obukhov predicted Du = 4 using an e-expansion 
renormalization group scheme pO| ]. Later Dfaz-Guilera 
performed a momentum space analysis of the correspond- 
ing Langevin equations which confirmed Du — 4 [ pd] |. 
Grassberger and Manna concluded from numerical in- 
vestigations of the BTW model in D < 5 the same result 
In contrast, comparable simulations and the simi- 
larity to percolation led several authors to the conjecture 
that Du = 6 comparable to the related forest fire 
model of Drossel and Schwabl (see |0 for an overview). 

In the present work we consider the BTW model in 
various dimensions {D < 6) on lattice sizes which are 
significant larger than those considered in previous works 
| p^ , p^JT5| . A finite size scaling analysis allows us to deter- 
mine the avalanche exponents, the dynamical exponent 
and to analyse whether the avalanche clusters are fractal. 
Our analysis reveals that the upper critical dimension is 



Du — 4 and that the avalanches display a fractal behav- 
ior above Du- We discuss the dimensional dependence 
of the exponents and derive scaling relations. Finally 
we briefly report results of similar investigations of the 
ZJ-state model which is a possible generalization of the 
two-state model introduced by Manna in two-dimensions 
|]l6| . It is known that the BTW model and Manna's 
model belong to different universality classes in D = 2 

Hi- 



ll. MODEL AND SIMULATIONS 

We consider the D-dimensional BTW model on a 
square lattice of linear size L in which integer variables 
hr > represent local heights. One perturbes the system 
by adding particles at a randomly chosen site accord- 
ing to 



1 



with random r. 



(1) 



A site is called unstable if the corresponding height 
exceeds a critical value he, i.e., if > he, where he is 
given by he = 2D. An unstable site relaxes, its value is 
decreased by he and the 2D next neighboring sites are 
increased by one unit, i.e., 

hr ^ hr — he (2) 
hnn,v ^ hnn,v ~t- 1. (3) 

In this way the neighboring sites may be activated and an 
avalanche of relaxation events may take place. The sites 
are updated in parallel until all sites are stable. Then 
the next particle is added [Eq. ([|)]. We assume open 
boundary conditions with heights at the boundary fixed 
to zero. 

System sizes L < 256 for D = 3, L < 80 for D = 4, 
L < 36 for D — b, and L < 18 for 13 = 6 are inves- 
tigated. Starting with a lattice of randomly distributed 
heights h e {0, 1,2, he—l} the system is perturbed ac- 
cording to Eq. (P and Dhar's 'burning algorithm' is ap- 
plied in order to check if the system has reached the crit- 
ical steady state [|]. Then we start the actual measure- 
ments which are averaged over at least 2 x lO'' non-zero 
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avalanches. We studied four different properties charac- 
terizing an avalanche: the number of relaxation events s, 
the number of distinct toppled lattice site sj. (area), the 
duration t, and the radius r. For a detailed description 
see |6| and references therein. In the critical steady state 
the corresponding probability distributions should obey 
power-law behavior characterized by exponents t^, r^, Tt 
and Tr according to 



Pr{r) - r-''-. 



(4) 
(5) 
(6) 
(7) 



Because a particular lattice site may topple several times 
the number of toppling events exceeds the number of dis- 
tinct toppled lattice sites, i.e., s > Sd- We will see that 
these multiple toppling events can be neglected for D > 3 
and the distribution Ps(s) and Pd{sd) display the same 
scaling behavior. 

Scaling relations for the exponents Ts,Td, Tt and Tr can 
be obtained if one assumes that the size, area, duration 
and radius scale as a power of each other, for instance 



t 



(8) 



The relation Pt{t)At — Pr{r)dr for the corresponding dis- 
tribution functions then leads to the scaling relation 



1 



Itr 



Tt-l 



(9) 



The exponents 7dr, 7rs, Isd etc are defined in the same 
way. The exponent is usually identified with the dy- 
namical exponent z and various theoretical efforts have 
been performed to determine z Dfaz-Guilera 
[ prf concluded from a momentum-space analysis of the 
corresponding Langevin equations that the dynamical ex- 
ponent of the BTW model is given by 



D 



(10) 



which was already suggested by Zhang Numerical 
investigations suggest that Eq. (^ is valid On 
the other hand Majumdar and Dhar used the equiv- 
alence between the sandpile model and the q — ^ limit 
of the Potts model to estimate z = | in _D = 2 which 
contradicts Eq. (p^). 

Christensen and Olami showed that inside an 
avalanche no holes can occur in the steady state jl^ ] 
where a hole is a set of untoppled sites which are com- 
pletely enclosed by toppled lattice sites. This implies for 
D — 2 that the avalanches are simply connected and com- 
pact. For D > 2 holes are still forbidden in the steady 
state but loops of toppled sites can occur. Then the 
avalanches are no more simply connected (see below). 



Even though no holes inside an avalanche cluster can oc- 
cur it was already assumed that above the critical dimen- 
sion Du the avalanches have the fractal dimension 4 
Here, the propagation of an avalanche can not be con- 
sidered as a connected activation front of toppled sites. 
The behavior is similar to an branching process where 
disconnected arms propagate without forming loops. If 
the avalanche clusters are not fractal the scaling exponent 
jdr which describes how the number of toppled sites Sd 
scales with the radius r equals the dimension D. Thus, 
the dimensional dependence of the exponent 7^,. is an 
appropriate tool to investigate the developing fractal be- 
havior with increasing dimension. 

The measurement of the probability distributions and 
the corresponding exponents [Eq. (Q^)] is affected by the 
finite systems size. For instance, the two dimensional 
BTW model displays a logarithmic wstem size depen- 
dence of the distribution exponents [jl^,^ . Another ex- 
ample is the related two dimensional Zhang model |l7| ] 
where the exponents depend on the inverse system size, 
i.e., the corrections are of the relative magnitude of the 
boundary L^^ In these cases the exponent of the 

infinite system could be obtained by an extrapolation to 
the infinite system size. If the values of the avalanche 
exponents r are not affected by the finite system size the 
powerful method of finite size scaling would be applica- 
ble. Here, the probability distributions [Eq. (^-0)] obey 
the scaling equation 



P,(x,L) = L-/^- g^iL-"- x), 



(11) 



with X £ js, d, t, r} and where is called the universal 
function pl| |. The exponent Tx is related to the scaling 
exponents f3x and Vx via 



(12) 



The exponent Vx determines the cut-off behavior of the 
probability distribution. If finite size scaling works all 
distributions Px{x,L) for various system sizes have to 
collapse, including their cut-offs. Then the argument 
of the universal function g has to be constant, i.e., 
XmaxL'"^ = const. Using the corresponding scahng re- 
lation [Eq. (||)] yields r^axL'"" = const. The cut-off 
radius r^ax should scale with the system size L and fi- 
nally one gets 



Ix 



(13) 



The advantage of the finite size scaling analysis is that 
it yields additionally to the avalanche exponents Tx the 
important scaling exponents ^dr and 74^ = z. 



III. D = 3 

In _D = 3 multiple toppling events, i.e., s > s^, oc- 
curs for less than 5% of all avalanches (nearly 42% in 
D ^ 2 and less than 0.1% in D = 4). These multiple 
toppling avalanches do not affect the scaling behavior of 
the probability distribution Ps{s)^ in the sense that there 
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is no significant difference between Ps{s) and Pd{sd) (see 
Fig. 1). Thus one concludes that tj^ = which is con- 
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FIG. 3. The scaling plot of the probability distribution 
Pd(sd) for L e {64,96, 128, 256} and D = 3. The dashed 
line corresponds to a power- law with an exponent |. 



FIG. 1. The probability distributions Ps{s) and Pd{sd) for 
L G {16, 32, 64, 128, 256}. For L < 256 the curves are shifted 
in the downward direction. Note that there is no significant 
difference between both distributions, i.e., multiple toppling 
events can be neglected in _D = 3. 
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FIG. 2. The system size dependence of the avalanches ex- 
ponents Td, Tt, and Tr for D = 3. The horizontal dashed lines 
correspond to the values |, |, and 2. The dotted lines are 
fits according to the equations Tx{L) = Tx ~ constL"^. 

The exponents t^, r^, and t^, obtained from a power- 
law fit of the straight portion of the probability distribu- 
tions [Eq. are plotted in Fig. g| for various system 
sizes L. The system size dependence vanishes quickly 
with increasing L. The dotted lines in Fig. ^ corresponds 
to a L^^ dependence of the avalanche exponents. The 
finite size corrections are of the magnitude of the bound- 
ary term in three dimensions. For L > 64 the system size 
dependence of Td and Tt is smaller than the statistical er- 
ror of the determination and the average of the exponents 
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FIG. 4. The scaling plot of the probability distribution 
Pt{t) for L G {64, 96, 128, ...,256} and D = 3. The dashed 
line corresponds to a power-law with an exponent | . 

for L > 64 would be a good estimate of the values of the 
infinite system. We obtain the values — 1.333 ± 0.007 
and Tt — 1.597 ± 0.012. The value of Td is in agree- 
ment with previous investigations based on smaller sys- 
tem sizes 1^,^. The exponent r,. seems to converge in 
the vicinity of 2 but the accuracy of this measurement is 
not sufficient to decide whether the value is exactly two. 
However, the following analysis lead us to the conclusion 
that Tr = 2. 

Since the avalanche exponents Td and Tt display no 
significant system size dependence for L > 64 the 
above mentioned finite-size scaling analysis is applicable 
[Eq. (pT|)]. The scaling plots of the distributions Pd{sd) 
and Ptit) are shown in Fig. ^ and Fig. IJ. One obtains 
a convincing data collapse of the various curves corre- 
sponding to the different system sizes for (3d = 4, = 3, 
and /3t = |, ft = z = |, respectively. Using Eq. (|l^ ) the 
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avalanches exponents are given by = | and Tt = |. 
These values are in agreement with our results obtained 
from a direct determination of the exponents via regres- 
sion. The value z = ^ agrees with Eq. (10) and Vd — 5 
reflects the fact that the avalanches are not fractal. This 
does not mean that the avalanche clusters are still sim- 
ply connected since the avalanches can form loops. But 
these rare loops do not contribute to the scaling behav- 
ior. Both scaling relations, (r^ — 1) = z{Tt — 1) and 
(r^ — 1) = jdr{Td — l), confirm our assumption that = 2. 
In summary our direct measurements as well as the finite 
size scaling analysis both yield that the avalanche expo- 
nents of the three dimensional BTW model are consistent 
with the values Td — Tg — |, — ^, Tr = 2, z — |, and 
Idr = 3. All scaling relations which connect these expo- 
nents are fulfilled. 



IV. D > 4 

Focusing our attention to the area and duration prob- 
ability distribution we find that finite size scaling works 
quite well again. In Fig. ^, ^, and ^ we present the scal- 
ing plots of the avalanche distribution Pd{sd) for D = A, 
D — 5, and D = 6. In all cases one gets a satisfying data 
collapse for Pd — G and Vd — 4, i.e., the corresponding 
avalanches exponent equals the mean field value = |. 
A similar analysis displays that the scaling exponents of 
the duration distribution Pt{t) are given by /3t = 4 and 
j/t = 2 resulting in Tt = 2 (not shown). The avalanche 
exponents of the BTW model in Z? > 4 agree with the 
mean field exponents Td = ^, Tt ~ 2, z = 2, and the up- 
per critical dimension is Z?„ = 4. All exponents are listed 
in Table | An analysis of the probability distribution 
Pr(r) and the determination of the exponent remains 
outside the scope of this paper because the considered 
system sizes (limited by computer power) are too small. 
For instance, in the case oi D = 4 the largest considered 
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FIG. 7. The scaling plot of the probability distribution 
Pdisd) for L e {8, 10, 12, 18} and D = 6. The dashed 
line corresponds to a power-law with the exponent | . 



system sizes is i = 80. The corresponding distribution 
Prir) exhibits a very small power-law region (less than 
a half decade), forbidding any accurate determination of 

Tr- 

The value jdr — 4 corresponds to the fact that the 
avalanches of the BTW model display a fractal behavior 
above the critical dimension Du , whereby the area scales 
with the radius according to Sd ^ r^, independently of 
the embedding dimension D. For D < Du the avalanches 
are not fractal. We display this developing fractal be- 
havior in Fig. 1^, where four arbritrary chosen avalanche 
clusters are shown for three different dimensions. For 
D > 4 we plotted three dimensional cuts through the 
center of mass of the avalanche clusters. The isolated is- 
lands which appear in the avalanche snapshots for D > 4 
are caused by the three dimensional cuts. In all cases 
the system size is L = 32 and the area of the plotted 
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FIG. 8. Snapshots of four arbritrary chosen avalanche clusters for D — 3 (upper left), D = 4 (upper right), and D — 5 
(lower left and right). For D > 4 three dimensional cuts through the center of mass are shown. In the three dimensional 
avalanche a loop can be seen in the upper left part of the avalanche cluster. 



avalanches is Sd = 1520 in _D = 3, Sd = 17500 in _D = 4, 
and Sd = 201000 in _D = 5, i.e., s^J''^ is nearly fixed. If 
the avalanches are not fractal in all dimensions the scal- 
ing relation Sd ~ y'' holds for all D and the radius should 
be independent of the embedding dimension. One can see 
see from Fig. ^ that the radius of the shown avalanches 
is roughly the same for D = 3 and D = A. Despite some 
loops (e.g. in the upper left part of the plotted three di- 
mensional avalanche) the avalanche clusters look nearly 
compact. In the five dimensional case the clusters dis- 
play a fractal behavior. The radius seems to be larger 
compared to the lower dimensional cases indicating that 
the equation Sd ~ does not hold in D = 5. Of course 
these snapshots only illustrate the developing fractal be- 
havior. 



Our results are in contrast to previous investigations 
performed by Janosi and Czirok |20|. They calculated 
the number of toppled site N{r) inside a sphere with ra- 
dius r. The sphere is centered at the center of mass of 
the avalanche cluster. The fractal dimension Df is ob- 



tained from the scaling law N{r) 



Considering one 



system size {L = 100 in Z? = 3) they found that the frac- 
tal dimension is given by -D/ « 2.75, i.e., the avalanches 
already display a fractal behavior in three dimensions. 
We performed the same analysis and reproduced their 
results within the error-bars. Analyzing various system 
sizes, however, we find that the apparent fractal dimen- 
sion depends on the system size and tends to Df = 3 with 
increasing L (not shown) in agreement with our results, 
discussed above. 
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TABLE I. Values of the exponents of the BTW model in 
various dimensions. 
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V. DISCUSSION 

In the following we examine the avalanche exponents 
as a function of the dimension D. Consider the average 
avalanche size 



y sPsis,L) ds. 



Using the finite size scaling ansatz [Eq. (|1 
for > 3 one gets 



(14) 



which works 



{s)i 



L 



(15) 



if Ts < 2. On the other hand it is known exactly g] 
that (s)l ^ in Z? = 2 and arguing that in undirected 
models particles diffuses out to the boundary one gets 
the same result independent of the dimension ||2l[] . Like 
Grassberger and Manna ||l2| we plot in Fig. ^ the average 
avalanches size as a function of the system size for var- 
ious dimensions. Except of deviations for small system 
sizes all data points collapse on a single curve. Thus one 
concludes that the equation 2 = 7sr(2 — t^) is fulfilled. 
Neglecting multiple toppling (ts = Td and 7sr = Idr) 
which is valid for D > i and using that the avalanches 
are not fractal {'-^dr — D) which is fulfilled for D < Du 
one gets 



Td = 2 , 



(16) 



for 3 < D < Du |2^]. This equation was already derived 
in the continuum limit by Zhang using energy conserva- 
tion and the local nature of energy transfer [n7| . Now we 
see that the failure of this equation for I? = 2 is caused 
by multiple toppling events which are essential in the two 
dimensional model only. For D > 3 multiple toppling can 
be neglected and Eq. (|6|) is fulfilled. Using 

z{Tt-l) = 7rfr(Td-l) (17) 

and Eq. ( px| ) the duration exponent Tt is given by 

n - 1 



FIG. 9. The average avalanche size (s) l as a function of the 
system size in various dimensions. The solid line corresponds 
to the power-law (s) l ^ which is exactly in D = 2 



again for 3 < D < Du- 

Finally we briefly report results of similar investiga- 
tions of the related ZJ-state sandpile model based on 
Manna's two-dimensional two-state model . Here the 
critical height he equals the dimension D and an unstable 
site relaxes to zero, whereby the particles are distributed 
randomly among the nearest neighbors. Again we find 
that the upper critical dimension is Du = 4. In contrast 
to the BTW model the dimensional dependence of the 
dynamical exponent is given by z = (Z? -t- 4)/4. Our pre- 
liminary results for D ~ "i are Ts ~ ^ , Td ~ ^ , Tt 
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10 ' 

Tr ~ |, and "fdr ~ 3. We find that Td is definitely larger 
than Ts (in agreement with |]l5|), i.e., multiple toppling 
events are relevant in the three dimensional model. Be- 
cause in the ZJ-state model the toppling processes are 
isotropic on average only holes inside an avalanche clus- 
ter can occur. But nevertheless, we find that ^dr = D for 
D < Du, i.e., these holes occur only on finite sizes and do 
not contribute to the scaling behavior. Above the critical 
dimension Du = 4 the avalanches have fractal dimension 
4. In D = 4 and D = 5 the model is characterized by the 
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mean field exponents, comparable to the BTW model. 
The values of the exponents are listed in Table ||. 

VI. CONCLUSIONS 

We studied numerically the dynamical properties of 
the BTW model on a square lattice in various dimen- 
sions. Using a finite size scaling analysis we determined 
the probability distribution exponents, the dynamical ex- 
ponent, and the dimension of the avalanches. Our anal- 
ysis reveals that multiple toppling events are relevant in 
the low dimensional case only and can be neglect for 
D > 3. For D = 3 the exponents are given by = 2, 
Tt = |, Trf = |, and 2 = §• For D > A the exponents 
agree with the mean field and Bethe lattice exponents, re- 
spectively. We conclude from our numerical results that 
below the critical dimension the dynamical exponent z 
is given by z = (D + 2)/3. The avalanche clusters are 
simply connected for £> = 2 only. For D > 2 loops occur 
but do not contribute to the scaling behavior until the 
embedding dimension exceeds the upper critical dimen- 
sion Du- Above Du the avalanches are fractal with the 
fractal dimension 4. 
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